#ifndef SERIESLEVIN_H
#define SERIESLEVIN_H

#include "SeriesBase.h"

// A series extrapolator using the Levin u transform

template <class T> 
class SeriesLevin : 
    public virtual SeriesBase<T>
{
#define INITVECS(n) \
    do { \
        q_num.reserve(n); \
        q_den.reserve(n); \
        if (this->doerrcalc) { \
            dq_num.reserve(n); \
            dq_den.reserve(n); \
            dsum.reserve(n); \
        } \
    } while (false)
public:
    SeriesLevin(uint n=20, bool _doerrcalc=true, bool _saveterms=true) : 
        SeriesBase<T>(n,_doerrcalc,_saveterms), result_n(0.), startsize(n)
    { INITVECS(n); }

    SeriesLevin(const std::vector<T>& v, bool _doerrcalc=true,
                bool _saveterms=true) : 
        SeriesBase<T>(v.size(),_doerrcalc,_saveterms), result_n(0.),
        startsize(v.size()) 
    { INITVECS(v.size()); AddTerms(v); }

    SeriesLevin(const std::vector<T>& v, const std::vector<T>& verr, 
                bool _doerrcalc=true, bool _saveterms=true) : 
        SeriesBase<T>(v.size(),_doerrcalc,_saveterms), result_n(0.),
        startsize(v.size()) 
    { INITVECS(v.size()); AddTerms(v,verr); }

    ~SeriesLevin() {}
#undef INITVECS

protected:
    virtual void CalcResult(const T& term);
    virtual const T& GetResult() const {return result_n;}
    virtual const std::vector<T>& GetDSum() const 
    {
        if (dsum.size() != q_num.size()) {
            uint n = dq_num[0].size();
            dsum.resize(n);
            assert(dq_den[0].size() == n);
            T result_n = GetResult();
            for (uint i = 0; i < n; i++) {
                dsum[i] = (dq_num[0][i] - result_n * dq_den[0][i]) / q_den[0];
            }
        }
        return dsum;
    }


private:
    std::vector<T> q_num; // backward diagonal of numerator
    std::vector<T> q_den; // backward diagonal of denominator
    std::vector<std::vector<T> > dq_num; // table of numerator derivatives
    // dq_num[i][j] = d(q_num[i])/d(term_j)
    std::vector<std::vector<T> > dq_den; // table of denominator derivatives
    // dq_den[i][j] = d(q_den[i])/d(term_j)
    mutable std::vector<T> dsum;
    T result_n;
    uint startsize;
};

template <class T>
inline void SeriesLevin<T>::CalcResult(const T& term)
{
    uint n = this->NTerms();
    assert(q_num.size()+1 == n);
    assert(q_den.size()+1 == n);
    if (this->doerrcalc) {
        assert(dq_den.size() == q_den.size());
        assert(dq_num.size() == q_den.size());
    }

    // Update the tables for the new term
    if (n == 1) {
        T invterm = T(1.0)/term;
        q_den.push_back(invterm);
        q_num.push_back(T(1.0));

        if (this->doerrcalc) {
            dq_den.push_back(std::vector<T>(1,-invterm*invterm));
            dq_den.back().reserve(startsize);
            dq_num.push_back(std::vector<T>(1,T(0.0)));
            dq_num.back().reserve(startsize);
        }
        result_n = term;
    } else {
        double nsqr = n*n;
        T invterm = 1.0/term;
        T sum_plain = this->PlainSum();
        q_den.push_back(invterm / nsqr);
        q_num.push_back(sum_plain * q_den[n-1]);

        if (this->doerrcalc) {
            for (uint i = 0; i < n-1; i++) {
                dq_den[i].push_back(T(0.0));
                dq_num[i].push_back(q_den[n-1]);
            }
            dq_den.push_back(std::vector<T>(n,-q_den[n-1]*invterm));
            dq_den.back().reserve(startsize);
            dq_num.push_back(
                std::vector<T>(n,q_den[n-1] + sum_plain*dq_den[n-1][n-1]));
            dq_num.back().reserve(startsize);
        }

        double factornum = 1.0;
        double factorden = n;

        for (int i = n-1; i > 0; i--) {
            T c = factornum * i;
            // If T is extra precision, this keeps all the precision with 
            // minimal calculation
            c /= factorden;
            q_den[i-1] = q_den[i] - c * q_den[i-1];
            q_num[i-1] = q_num[i] - c * q_num[i-1];

            if (this->doerrcalc) {
                for (uint j = 0; j < n-1; j++) {
                    dq_den[i-1][j] = dq_den[i][j] - c * dq_den[i-1][j];
                    dq_num[i-1][j] = dq_num[i][j] - c * dq_num[i-1][j];
                }
            }
            factornum *= double(n-1);
            factorden *= double(n);
        }
        result_n = q_num[0]/q_den[0];
    }
}

// Specialize for complex<T> Series.  Only difference is that some of 
// that variables for the calculations in AddTerm don't have to be 
// complex, so they are a little faster this way.

template <class T> 
class SeriesLevin<std::complex<T> > : 
    public virtual SeriesBase<std::complex<T> >
{
#define INITVECS(n) \
    do { \
        q_num.reserve(n); \
        q_den.reserve(n); \
        if (this->doerrcalc) { \
            dq_num.reserve(n); \
            dq_den.reserve(n); \
            dsum.reserve(n); \
        } \
    } while (false)
public:
    SeriesLevin(uint n=20, bool _doerrcalc=true, bool _saveterms=true) : 
        SeriesBase<std::complex<T> >(n,_doerrcalc,_saveterms),
        result_n(0.), startsize(n) 
    { INITVECS(n); }

    SeriesLevin(const std::vector<std::complex<T> >& v,
                bool _doerrcalc=true, bool _saveterms=true) : 
        SeriesBase<std::complex<T> >(v,_doerrcalc,_saveterms), result_n(0.), 
        startsize(v.size()) 
    { INITVECS(v.size()); }

    SeriesLevin(const std::vector<std::complex<T> >& v,
                const std::vector<std::complex<T> >& verr, 
                bool _doerrcalc=true, bool _saveterms=true) : 
        SeriesBase<std::complex<T> >(v,verr,_doerrcalc,_saveterms),
        result_n(0.), startsize(v.size())
    { INITVECS(v.size()); }

    ~SeriesLevin() {}
#undef INITVECS

protected:
    virtual void CalcResult(const std::complex<T>& term);
    virtual const std::complex<T>& GetResult() const {return result_n;}
    virtual const std::vector<std::complex<T> >& GetDSum() const 
    {
        if (dsum.size() != q_num.size()) {
            uint n = dq_num[0].size();
            dsum.resize(n);
            assert(dq_den[0].size() == n);
            std::complex<T>  result_n = GetResult();
            for (uint i = 0; i < n; i++) {
                dsum[i] = (dq_num[0][i] - result_n * dq_den[0][i]) / q_den[0];
            }
        }
        return dsum;
    }


private:
    std::vector<std::complex<T> > q_num; 
    std::vector<std::complex<T> > q_den; 
    std::vector<std::vector<std::complex<T> > > dq_num; 
    std::vector<std::vector<std::complex<T> > > dq_den;
    mutable std::vector<std::complex<T> > dsum;
    std::complex<T> result_n;
    uint startsize;
};

template <class T>
inline void SeriesLevin<std::complex<T> >::CalcResult(
    const std::complex<T>& term)
{
    typedef std::complex<T> CT;
    uint n = this->NTerms();
    assert(q_num.size()+1 == n);
    assert(q_den.size()+1 == n);
    if (this->doerrcalc) {
        assert(dq_den.size() == q_den.size());
        assert(dq_num.size() == q_den.size());
    }

    if (n == 1) {
        CT invterm = 1.0/term;
        q_den.push_back(invterm);
        q_num.push_back(CT(1.0));

        if (this->doerrcalc) {
            dq_den.push_back(std::vector<CT>(1,-invterm*invterm));
            dq_den.back().reserve(startsize);
            dq_num.push_back(std::vector<CT>(1,CT(0.0)));
            dq_num.back().reserve(startsize);
        }
        result_n = term;
    } else {
        double nsqr = n*n;
        CT invterm = 1.0/term;
        CT sum_plain = this->PlainSum();
        q_den.push_back(invterm / nsqr);
        q_num.push_back(sum_plain * q_den[n-1]);

        if (this->doerrcalc) {
            for (uint i = 0; i < n-1; i++) {
                dq_den[i].push_back(CT(0.0));
                dq_num[i].push_back(q_den[n-1]);
            }
            dq_den.push_back(std::vector<CT>(n,-q_den[n-1]*invterm));
            dq_den.back().reserve(startsize);
            dq_num.push_back(
                std::vector<CT>(n,q_den[n-1] + sum_plain*dq_den[n-1][n-1]));
            dq_num.back().reserve(startsize);
        }

        double factornum = 1.0;
        double factorden = n;

        for (int i = n-1; i > 0; i--) {
            T c = factornum * i;
            c /= factorden;
            q_den[i-1] = q_den[i] - c * q_den[i-1];
            q_num[i-1] = q_num[i] - c * q_num[i-1];

            if (this->doerrcalc) {
                for (uint j = 0; j < n-1; j++) {
                    dq_den[i-1][j] = dq_den[i][j] - c * dq_den[i-1][j];
                    dq_num[i-1][j] = dq_num[i][j] - c * dq_num[i-1][j];
                }
            }
            factornum *= double(n-1);
            factorden *= double(n);
        }
        result_n = q_num[0]/q_den[0];
    }
}

#endif
